Introduction à la modélisation statistique bayésienne

Planning

Cours n°01 : Introduction à l’inférence bayésienne
Cours n°02 : Modèle Beta-Binomial
Cours n°03 : Introduction à brms, modèle de régression linéaire
Cours n°04 : Modèle de régression linéaire (suite)
Cours n°05 : Markov Chain Monte Carlo
Cours n°06 : Modèle linéaire généralisé
Cours n°07 : Comparaison de modèles
Cours n°08 : Modèles multi-niveaux
Cours n°09 : Modèles multi-niveaux généralisés
Cours n°10 : Data Hackathon

Langage de la modélisation

\[ \begin{align} y_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i}&= \alpha + \beta x_{i} \\ \alpha &\sim \mathrm{Normal}(60, 10) \\ \beta &\sim \mathrm{Normal}(0, 10) \\ \sigma &\sim \mathrm{HalfCauchy}(0, 1) \end{align} \]

Objectif de la séance : comprendre ce type de modèle.

Un premier modèle

library(tidyverse)
library(imsb)

d <- get(data(howell) )
str(d)

'data.frame':   544 obs. of  4 variables:
 $ height: num  152 140 137 157 145 ...
 $ weight: num  47.8 36.5 31.9 53 41.3 ...
 $ age   : num  63 63 65 41 51 35 32 27 19 54 ...
 $ male  : int  1 0 0 1 0 1 0 1 0 1 ...

d2 <- d %>% filter(age >= 18)
head(d2)

   height   weight age male
1 151.765 47.82561  63    1
2 139.700 36.48581  63    0
3 136.525 31.86484  65    0
4 156.845 53.04191  41    1
5 145.415 41.27687  51    0
6 163.830 62.99259  35    1

Un premier modèle

\[h_{i} \sim \mathrm{Normal}(\mu, \sigma)\]

d2 %>%
    ggplot(aes(x = height) ) +
    geom_histogram(bins = 10, col = "white")

Loi normale

\[ p(x \ | \ \mu, \sigma) = \frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \bigg[-\frac{1}{2 \sigma^{2}} (\mu - x)^{2} \bigg] \]

data.frame(value = rnorm(1e4, 10, 1) ) %>% # 10.000 samples from Normal(10, 1)
    ggplot(aes(x = value) ) +
    geom_histogram(col = "white")

D’où vient la loi normale ?

Certaines valeurs sont fortement probables (autour de la moyenne \(\mu\)). Plus on s’éloigne, moins les valeurs sont probables (en suivant une décroissance exponentielle).

D’où vient la loi normale ?

\[ y = \exp \big[-x^{2} \big] \]

On étend notre fonction aux valeurs négatives.

D’où vient la loi normale ?

\[ y = \exp \big[-x^{2} \big] \]

Les points d’inflection nous donnent une bonne indication de là où la plupart des valeurs se trouvent (i.e., entre les points d’inflection). Les pics de la dérivée nous montrent les points d’inflection.

D’où vient la loi normale ?

\[ y = \exp \bigg [- \frac{1}{2} x^{2} \bigg] \]

Ensuite on standardise la distribution de manière à ce que les deux points d’inflection se trouvent à \(x = -1\) et \(x = 1\).

D’où vient la loi normale ?

\[ y = \exp \bigg [- \frac{1}{2 \color{steelblue}{\sigma^{2}}} x^{2} \bigg] \]

On insère un paramètre \(\sigma^{2}\) pour contrôler la distance entre les points d’inflection.

D’où vient la loi normale ?

\[ y = \exp \bigg [- \frac{1}{2 \color{steelblue}{\sigma^{2}}} (\color{orangered}{\mu} - x)^{2} \bigg] \]

On insère ensuite un paramètre \(\mu\) afin de pouvoir contrôler la position (la tendance centrale) de la distribution.

D’où vient la loi normale ?

\[ y = \frac{1}{\sqrt{2 \pi \color{steelblue}{\sigma^{2}}}} \exp \bigg[-\frac{1}{2 \color{steelblue}{\sigma^{2}}} (\color{orangered}{\mu} - x)^{2} \bigg] \]

Mais… cette distribution n’intègre pas à 1. On divise donc par une constante de normalisation (la partie gauche), afin d’obtenir une distribution de probabilité.

Modèle gaussien

Nous allons construire un modèle de régression, mais avant d’ajouter un prédicteur, essayons de modéliser la distribution des tailles.

Modèle gaussien

On définit \(p(\mu, \sigma)\), la distribution a priori conjointe de tous les paramètres du modèle. On peut spécifier ces priors indépendamment pour chaque paramètre, sachant que \(p(\mu, \sigma) = p(\mu) p(\sigma)\).

\[\color{steelblue}{\mu \sim \mathrm{Normal}(178, 20)}\]

Modèle gaussien

On définit \(p(\mu, \sigma)\), la distribution a priori conjointe de tous les paramètres du modèle. On peut spécifier ces priors indépendamment pour chaque paramètre, sachant que \(p(\mu, \sigma) = p(\mu) p(\sigma)\).

\[\color{steelblue}{\sigma \sim \mathrm{Uniform}(0, 50)}\]

Visualiser le prior

library(ks)
sample_mu <- rnorm(1e4, 178, 20) # prior on mu
sample_sigma <- runif(1e4, 0, 50) # prior on sigma
prior <- data.frame(cbind(sample_mu, sample_sigma) ) # multivariate prior
H.scv <- Hscv(x = prior, verbose = TRUE)
fhat_prior <- kde(x = prior, H = H.scv, compute.cont = TRUE)
plot(
    fhat_prior, display = "persp", col = "steelblue", border = NA,
    xlab = "\nmu", ylab = "\nsigma", zlab = "\n\np(mu, sigma)",
    shade = 0.8, phi = 30, ticktype = "detailed",
    cex.lab = 1.2, family = "Helvetica")

Prior predictive checking

sample_mu <- rnorm(1000, 178, 20)
sample_sigma <- runif(1000, 0, 50)

data.frame(x = rnorm(1000, sample_mu, sample_sigma) ) %>%
    ggplot(aes(x) ) +
    geom_histogram() +
    labs(x = "Taille (en cm)", y = "Nombre d'échantillons")

Fonction de vraisemblance

mu_exemple <- 151.23
sigma_exemple <- 23.42

d2$height[34] # une observation de taille (pour exemple)

[1] 162.8648

Fonction de vraisemblance

On veut calculer la probabilité d’observer une certaine valeur de taille, sachant certaines valeurs de \(\mu\) et \(\sigma\), c’est à dire :

\[ p(x \ | \ \mu, \sigma) = \frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \bigg[-\frac{1}{2 \sigma^{2}} (\mu - x)^{2} \bigg] \]

dnorm(d2$height[34], mu_exemple, sigma_exemple)

[1] 0.01505675

Fonction de vraisemblance

\[ p(x \ | \ \mu, \sigma) = \frac{1}{\sqrt{2 \pi \sigma^{2}}} \exp \bigg[-\frac{1}{2 \sigma^{2}} (\mu - x)^{2} \bigg] \]

Ou à la main…

normal_likelihood <- function (x, mu, sigma) {
  
  bell <- exp( (- 1 / (2 * sigma^2) ) * (mu - x)^2 )
  norm <- sqrt(2 * pi * sigma^2)
  
  return(bell / norm)
  
}

normal_likelihood(d2$height[34], mu_exemple, sigma_exemple)

[1] 0.01505675

Distribution postérieure

\[ \color{purple}{p(\mu, \sigma \ | \ h)} = \frac{\prod_{i} \color{orangered}{\mathrm{Normal}(h_{i} \ | \ \mu, \sigma)}\color{steelblue}{\mathrm{Normal}(\mu \ | \ 178, 20)\mathrm{Uniform}(\sigma \ | \ 0, 50)}} {\color{green}{\int \int \prod_{i} \mathrm{Normal}(h_{i} \ | \ \mu, \sigma)\mathrm{Normal}(\mu \ | \ 178, 20)\mathrm{Uniform}(\sigma \ | \ 0, 50) \mathrm{d} \mu \mathrm{d} \sigma}} \]

\[ \color{purple}{p(\mu, \sigma \ | \ h)} \propto \prod_{i} \color{orangered}{\mathrm{Normal}(h_{i} \ | \ \mu, \sigma)}\color{steelblue}{\mathrm{Normal}(\mu \ | \ 178, 20)\mathrm{Uniform}(\sigma \ | \ 0, 50)} \]

Distribution postérieure - Grid approximation

# définit une grille de valeurs possibles pour mu et sigma
mu.list <- seq(from = 140, to = 160, length.out = 200)
sigma.list <- seq(from = 4, to = 9, length.out = 200)

# étend la grille en deux dimensions (chaque combinaison de mu et sigma)
post <- expand.grid(mu = mu.list, sigma = sigma.list)

# calcul de la log-vraisemblance (pour chaque couple de mu et sigma)
post$LL <-
  sapply(
    1:nrow(post),
    function(i) sum(dnorm(
      d2$height,
      mean = post$mu[i],
      sd = post$sigma[i],
      log = TRUE)
      )
    )

# calcul de la probabilité a posteriori (non normalisée)
post$prod <-
  post$LL +
  dnorm(post$mu, 178, 20, log = TRUE) +
  dunif(post$sigma, 0, 50, log = TRUE)

# on "annule" le log en avec exp() et on standardise par la valeur maximale
post$prob <- exp(post$prod - max(post$prod) )

Distribution postérieure - Grid approximation

# select random 20 rows of the dataframe 
post %>% slice_sample(n = 20, replace = FALSE)

         mu    sigma        LL      prod          prob
1  147.0352 8.547739 -1360.471 -1369.496  3.269714e-62
2  154.2714 6.562814 -1230.413 -1238.943  1.632201e-05
3  156.4824 5.909548 -1267.972 -1276.377  9.026446e-22
4  154.5729 5.381910 -1279.103 -1287.616  1.187712e-26
5  141.1055 4.929648 -2636.175 -2645.703  0.000000e+00
6  157.1859 6.512563 -1258.869 -1267.237  8.415073e-18
7  149.3467 7.467337 -1306.848 -1315.701  7.540861e-39
8  159.4975 7.768844 -1289.440 -1297.694  4.983936e-31
9  148.6432 5.281407 -1510.091 -1518.995 3.871159e-127
10 158.5930 6.713568 -1289.476 -1297.774  4.603459e-31
11 142.8141 4.326633 -2706.399 -2715.773  0.000000e+00
12 157.7889 6.135678 -1289.110 -1297.448  6.377633e-31
13 158.9950 8.723618 -1268.879 -1277.157  4.137099e-22
14 148.5427 8.949749 -1306.803 -1315.714  7.443406e-39
15 156.4824 4.854271 -1352.575 -1360.980  1.632725e-58
16 148.2412 8.170854 -1326.938 -1335.872  1.310497e-47
17 150.7538 5.809045 -1331.574 -1340.329  1.519004e-49
18 141.3065 4.276382 -3110.214 -3119.724  0.000000e+00
19 147.4372 8.472362 -1347.877 -1356.871  9.938309e-57
20 148.1407 7.517588 -1349.505 -1358.446  2.057730e-57

Distribution postérieure - Grid approximation

sample.rows <- sample(x = 1:nrow(post), size = 1e4, replace = TRUE, prob = post$prob)

Distribution postérieure - Distributions marginales

posterior_plot(samples = sample.mu, nbins = 40) + labs(x = expression(mu) )

Distribution postérieure - Distributions marginales

posterior_plot(samples = sample.sigma, nbins = 40) + labs(x = expression(sigma) )

Introduction à brms

Under the hood : Stan est un langage de programmation probabiliste écrit en C++, et qui implémente plusieurs algorithmes de MCMC : HMC, NUTS, L-BFGS…

data {
  int<lower=0> J; // number of schools 
  real y[J]; // estimated treatment effects
  real<lower=0> sigma[J]; // s.e. of effect estimates 
  }

parameters {
  real mu; 
  real<lower=0> tau;
  real eta[J];
  }

transformed parameters {
  real theta[J];
  for (j in 1:J)
    theta[j] = mu + tau * eta[j];
  }

model {
  target += normal_lpdf(eta | 0, 1);
  target += normal_lpdf(y | theta, sigma);
  }

Bayesian regression models using Stan

Le package brms (Bürkner, 2017) permet de fitter des modèles multi-niveaux (ou pas) linéaires (ou pas) bayésiens en Stan mais en utilisant la syntaxe de lme4.

model <- brm(y ~ x + (1 | subject) + (1 | item), data = d, family = gaussian() )

Rappels de syntaxe

Le package brms utilise la même syntaxe que les fonctions de base R (comme lm) ou que le package lme4.

Reaction ~ Days + (1 + Days | Subject)

mvbind(Reaction, Memory) ~ Days + (1 + Days | Subject)

mvbind(Reaction, Memory) ~ Days + (1 + Days | Subject)
mvbind(Reaction, Memory) ~ 1 + Days + (1 + Days | Subject)

Rappels de syntaxe

Si l’on veut fitter un modèle sans intercept (why not), il faut le spécifier explicitement comme ci-dessous.

mvbind(Reaction, Memory) ~ 0 + Days + (1 + Days | Subject)

brm(Reaction ~ 1 + Days + (1 + Days | Subject), family = lognormal() )

Quelques fonctions utiles

# générer le code du modèle en Stan
make_stancode(formula, ...)
stancode(fit)

# définir les priors
get_prior(formula, ...)
set_prior(prior, ...)

# récupérer les prédiction du modèle
fitted(fit, ...)
predict(fit, ...)
conditional_effects(fit, ...)

# posterior predictive checking
pp_check(fit, ...)

# comparaison de modèles
loo(fit1, fit2, ...)
bayes_factor(fit1, fit2, ...)
model_weights(fit1, fit2, ...)

# test d'hypothèse
hypothesis(fit, hypothesis, ...)

Un premier exemple

library(brms)
mod1 <- brm(height ~ 1, data = d2)

posterior_summary(mod1, pars = c("^b_", "sigma"), probs = c(0.025, 0.975) )

              Estimate Est.Error       Q2.5      Q97.5
b_Intercept 154.599561 0.4099594 153.802456 155.409710
sigma         7.746558 0.2948896   7.208381   8.341852

En utilisant notre prior

Par défaut brms utilise un prior très peu informatif centré sur la valeur moyenne de la variable mesurée. On peut donc affiner l’estimation réalisée par ce modèle en utilisant nos connaissances sur la distribution habituelle des tailles chez les humains.

get_prior(height ~ 1, data = d2)

                    prior     class coef group resp dpar nlpar lb ub  source
 student_t(3, 154.3, 8.5) Intercept                                  default
     student_t(3, 0, 8.5)     sigma                             0    default

En utilisant notre prior

priors <- c(
  prior(normal(178, 20), class = Intercept),
  prior(exponential(0.01), class = sigma)
  )

mod2 <- brm(
  height ~ 1,
  prior = priors,
  family = gaussian(),
  data = d2
  )

En utilisant notre prior

summary(mod2)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ 1 
   Data: d2 (Number of observations: 352) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   154.61      0.42   153.78   155.45 1.00     3727     2790

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     7.78      0.30     7.21     8.39 1.00     2935     2564

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

En utilisant un prior plus informatif

priors <- c(
  prior(normal(178, 0.1), class = Intercept),
  prior(exponential(0.01), class = sigma)
  )

mod3 <- brm(
  height ~ 1,
  prior = priors,
  family = gaussian(),
  data = d2
  )

En utilisant un prior plus informatif

summary(mod3)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ 1 
   Data: d2 (Number of observations: 352) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept   177.86      0.10   177.66   178.06 1.00     3340     2404

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma    24.61      0.94    22.88    26.50 1.00     3666     2859

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

On remarque que la valeur estimée pour \(\mu\) n’a presque pas “bougée” du prior…mais on remarque également que la valeur estimée pour \(\sigma\) a largement augmentée. Nous avons dit au modèle que nous étions assez certain de notre valeur de \(\mu\), le modèle s’est ensuite “adapté”, ce qui explique la valeur de \(\sigma\)…

Précision du prior (heuristique)

Le prior peut généralement être considéré comme un posterior obtenu sur des données antérieures.

Visualiser les échantillons de la distribution postérieure

post <- as_draws_df(x = mod2) %>%
    mutate(density = get_density(b_Intercept, sigma, n = 1e2) )

ggplot(post, aes(x = b_Intercept, y = sigma, color = density) ) +
    geom_point(size = 2, alpha = 0.5, show.legend = FALSE) +
    labs(x = expression(mu), y = expression(sigma) ) +
    viridis::scale_color_viridis()

Récupérer les échantillons de la distribution postérieure

# gets the first 6 samples
head(post)

# A draws_df: 6 iterations, 1 chains, and 5 variables
  b_Intercept sigma lprior  lp__ density
1         154   7.4   -9.3 -1228   0.200
2         155   8.1   -9.3 -1228   0.312
3         155   8.0   -9.3 -1227   0.560
4         154   8.5   -9.3 -1231   0.011
5         154   7.4   -9.3 -1228   0.218
6         154   7.4   -9.3 -1228   0.386
# ... hidden reserved variables {'.chain', '.iteration', '.draw'}

# gets the median and the 95% credible interval
t(sapply(post[, 1:2], quantile, probs = c(0.025, 0.5, 0.975) ) )

                  2.5%        50%      97.5%
b_Intercept 153.775820 154.600133 155.447522
sigma         7.208968   7.765338   8.386128

Visualiser la distribution postérieure

H.scv <- Hscv(post[, 1:2])
fhat_post <- kde(x = post[, 1:2], H = H.scv, compute.cont = TRUE)

plot(fhat_post, display = "persp", col = "purple", border = NA,
  xlab = "\nmu", ylab = "\nsigma", zlab = "\np(mu, sigma)",
  shade = 0.8, phi = 30, ticktype = "detailed",
  cex.lab = 1.2, family = "Helvetica")

Visualiser la distribution postérieure

Ajouter un prédicteur

Comment est-ce que la taille co-varie avec le poids ?

d2 %>%
  ggplot(aes(x = weight, y = height) ) +
  geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8)

Régression linéaire à un prédicteur continu

\[ \begin{align} h_{i} &\sim \mathrm{Normal}(\mu_{i}, \sigma) \\ \mu_{i} &= \alpha + \beta x_{i} \\ \end{align} \]

linear_model <- lm(height ~ weight, data = d2)
rethinking::precis(linear_model, prob = 0.95)

                   mean         sd        2.5%       97.5%
(Intercept) 113.8793936 1.91106523 110.1337746 117.6250126
weight        0.9050291 0.04204752   0.8226175   0.9874407

Différentes notations

On considère un modèle de régression linéaire avec un seul prédicteur, une pente, un intercept, et des résidus distribués selon une loi normale. La notation :

\[ h_{i} = \alpha + \beta x_{i} + \epsilon_{i} \quad \text{avec} \quad \epsilon_{i} \sim \mathrm{Normal}(0, \sigma) \]

est équivalente à :

\[ h_{i} - (\alpha + \beta x_{i}) \sim \mathrm{Normal}(0, \sigma) \]

Régression linéaire à un prédicteur continu

\[ \begin{aligned} \color{orangered}{h_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i},\sigma)} \\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha + \beta x_{i}} \\ \color{steelblue}{\alpha} \ &\color{steelblue}{\sim \mathrm{Normal}(178, 20)} \\ \color{steelblue}{\beta} \ &\color{steelblue}{\sim \mathrm{Normal}(0, 10)} \\ \color{steelblue}{\sigma} \ &\color{steelblue}{\sim \mathrm{Exponential}(0.01)} \\ \end{aligned} \]

Dans ce modèle \(\mu\) n’est plus un paramètre à estimer (car \(\mu\) est déterminé par \(\alpha\) et \(\beta\)). À la place, nous allons estimer \(\alpha\) et \(\beta\).

Régression linéaire à un prédicteur continu

priors <- c(
  prior(normal(178, 20), class = Intercept),
  prior(normal(0, 10), class = b),
  prior(exponential(0.01), class = sigma)
  )

mod4 <- brm(
  height ~ 1 + weight,
  prior = priors,
  family = gaussian(),
  data = d2
  )

Régression linéaire à un prédicteur continu

posterior_summary(mod4)

                 Estimate  Est.Error          Q2.5         Q97.5
b_Intercept   113.8757806 1.93667879   110.1033438   117.7354033
b_weight        0.9052745 0.04265862     0.8198219     0.9885244
sigma           5.1067011 0.19087059     4.7491627     5.4997448
lprior        -12.4808134 0.01577149   -12.5112057   -12.4489637
lp__        -1083.3551723 1.24083628 -1086.6701725 -1081.9635005

• \(\alpha = 113.87, 95\% \ \text{CrI} \ [110.18, 117.68]\) représente la taille moyenne quand le poids est égal à 0kg…

• \(\beta = 0.91, 95\% \ \text{CrI} \ [0.82, 0.99]\) nous indique qu’une augmentation de 1kg entraîne une augmentation de 0.90cm.

Régression linéaire à un prédicteur continu

d2$weight.c <- d2$weight - mean(d2$weight)

mod5 <- brm(
  height ~ 1 + weight.c,
  prior = priors,
  family = gaussian(),
  data = d2
  )

fixef(mod5) # retrieves the fixed effects estimates

             Estimate  Est.Error        Q2.5       Q97.5
Intercept 154.6075906 0.27562525 154.0771202 155.1603670
weight.c    0.9049348 0.04264384   0.8231796   0.9886077

Représenter les prédictions du modèle

d2 %>%
    ggplot(aes(x = weight, y = height) ) +
    geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
    geom_abline(intercept = fixef(mod4)[1], slope = fixef(mod4)[2], lwd = 1)

Représenter l’incertitude sur \(\mu\) via fitted()

# on crée un vecteur de valeurs possibles pour "weight"
weight.seq <- data.frame(weight = seq(from = 25, to = 70, by = 1) )

# on récupère les prédictions du modèle pour ces valeurs de poids
mu <- data.frame(fitted(mod4, newdata = weight.seq) ) %>% bind_cols(weight.seq)

# on affiche les 10 premières lignes de mu
head(mu, 10)

   Estimate Est.Error     Q2.5    Q97.5 weight
1  136.5076 0.8926701 134.7786 138.2705     25
2  137.4129 0.8520675 135.7571 139.0942     26
3  138.3182 0.8116758 136.7355 139.9109     27
4  139.2235 0.7715282 137.7127 140.7401     28
5  140.1287 0.7316649 138.6915 141.5625     29
6  141.0340 0.6921348 139.6731 142.3849     30
7  141.9393 0.6529987 140.6532 143.2095     31
8  142.8446 0.6143317 141.6278 144.0336     32
9  143.7498 0.5762284 142.6105 144.8609     33
10 144.6551 0.5388083 143.5887 145.6965     34

Représenter l’incertitude sur \(\mu\) via fitted()

d2 %>%
  ggplot(aes(x = weight, y = height) ) +
  geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
  geom_smooth(
    data = mu, aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
    stat = "identity",
    color = "black", alpha = 0.8, size = 1
    )

Intervalles de prédiction (incorporer \(\sigma\))

Pour rappel, voici notre modèle : \(h_{i} \sim \mathrm{Normal}(\alpha + \beta x_{i}, \sigma)\). Pour l’instant, on a seulement représenté les prédictions pour \(\mu\). Comment incorporer \(\sigma\) dans nos prédictions ?

# on crée un vecteur de valeurs possibles pour "weight"
weight.seq <- data.frame(weight = seq(from = 25, to = 70, by = 1) )

# on récupère les prédictions du modèle pour ces valeurs de poids
pred_height <- data.frame(predict(mod4, newdata = weight.seq) ) %>% bind_cols(weight.seq)

# on affiche les 10 premières lignes de pred_height
head(pred_height, 10)

   Estimate Est.Error     Q2.5    Q97.5 weight
1  136.4999  5.179251 126.2393 146.6784     25
2  137.3616  5.088772 127.4212 147.8064     26
3  138.2271  5.156854 128.0879 148.1950     27
4  139.1509  5.177190 129.1437 149.6309     28
5  140.1046  5.050287 130.2454 149.8297     29
6  141.0736  5.207840 130.9432 151.4189     30
7  141.9294  5.090422 131.7089 151.8672     31
8  142.7915  5.104629 132.8222 152.8109     32
9  143.9081  5.181164 133.6667 154.0952     33
10 144.7201  5.155773 134.6358 155.0539     34

Intervalles de prédiction (incorporer \(\sigma\))

d2 %>%
  ggplot(aes(x = weight, y = height) ) +
  geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
  geom_ribbon(
    data = pred_height, aes(x = weight, ymin = Q2.5, ymax = Q97.5),
    alpha = 0.2, inherit.aes = FALSE
    ) +
  geom_smooth(
    data = mu, aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
    stat = "identity", color = "black", alpha = 0.8, size = 1
    )

Deux types d’incertitude

Deux sources d’incertitude dans le modèle : incertitude concernant l’estimation de la valeur des paramètres mais également concernant le processus d’échantillonnage.

Régression polynomiale

d %>% # on utilise d au lieu de d2
  ggplot(aes(x = weight, y = height) ) +
  geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8)

Scores standardisés

d <- d %>% mutate(weight.s = (weight - mean(weight) ) / sd(weight) )

d %>%
    ggplot(aes(x = weight.s, y = height) ) +
    geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8)

c(mean(d$weight.s), sd(d$weight.s) )

[1] -2.712698e-18  1.000000e+00

Scores standardisés

Pourquoi standardiser les prédicteurs ?

• Interprétation. Permet de comparer les coefficients de plusieurs prédicteurs. Un changement d’un écart-type du prédicteur correspond à un changement d’un écart-type sur la réponse (si la réponse est aussi standardisée).

• Fitting. Quand les prédicteurs contiennent de grandes valeurs, cela peut poser des problèmes de convergence (cf. Cours n°05)…

Modèle de régression polynomiale - exercice

\[ \begin{aligned} \color{orangered}{h_{i}} \ &\color{orangered}{\sim \mathrm{Normal}(\mu_{i}, \sigma)} \\ \color{black}{\mu_{i}} \ &\color{black}{= \alpha + \beta_{1} x_{i} + \beta_{2} x_{i}^{2}} \\ \color{steelblue}{\alpha} \ &\color{steelblue}{\sim \mathrm{Normal}(156, 100)} \\ \color{steelblue}{\beta_{1}, \beta_{2}} \ &\color{steelblue}{\sim \mathrm{Normal}(0, 10)} \\ \color{steelblue}{\sigma} \ &\color{steelblue}{\sim \mathrm{Exponential}(0.01)} \\ \end{aligned} \]

À vous de construire et fitter ce modèle en utilisant brms::brm().

Modèle de régression polynomiale

priors <- c(
  prior(normal(156, 100), class = Intercept),
  prior(normal(0, 10), class = b),
  prior(exponential(0.01), class = sigma)
  )

mod6 <- brm(
  # NB: polynomials should be written with the I() function...
  height ~ 1 + weight.s + I(weight.s^2),
  prior = priors,
  family = gaussian(),
  data = d
  )

Modèle de régression polynomiale

summary(mod6)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ 1 + weight.s + I(weight.s^2) 
   Data: d (Number of observations: 544) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     146.66      0.37   145.91   147.39 1.00     3466     2797
weight.s       21.40      0.29    20.83    21.96 1.00     3524     3016
Iweight.sE2    -8.41      0.28    -8.97    -7.86 1.00     3699     2893

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     5.78      0.18     5.44     6.14 1.00     3709     2854

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Représenter les prédictions du modèle

# on crée un vecteur de valeurs possibles pour "weight"
weight.seq <- data.frame(weight.s = seq(from = -2.5, to = 2.5, length.out = 50) )

# on récupère les prédictions du modèle pour ces valeurs de poids
mu <- data.frame(fitted(mod6, newdata = weight.seq) ) %>% bind_cols(weight.seq)
pred_height <- data.frame(predict(mod6, newdata = weight.seq) ) %>% bind_cols(weight.seq)

# on affiche les 10 premières lignes de pred_height
head(pred_height, 10)

   Estimate Est.Error     Q2.5     Q97.5  weight.s
1  40.80996  5.940786 29.42279  52.56775 -2.500000
2  46.84551  6.014583 35.34829  58.55580 -2.397959
3  53.39215  5.757048 42.16915  64.82299 -2.295918
4  59.33112  5.869666 48.05998  70.92279 -2.193878
5  65.12877  5.706725 53.66004  76.26556 -2.091837
6  70.80272  5.788879 59.54600  82.22328 -1.989796
7  76.39909  5.861752 64.87660  87.95047 -1.887755
8  81.71582  5.755260 70.50958  93.00193 -1.785714
9  86.81933  5.806592 75.37122  98.01592 -1.683673
10 91.87707  5.820529 80.33321 103.22341 -1.581633

Représenter les prédictions du modèle

d %>%
  ggplot(aes(x = weight.s, y = height) ) +
  geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
  geom_ribbon(
    data = pred_height, aes(x = weight.s, ymin = Q2.5, ymax = Q97.5),
    alpha = 0.2, inherit.aes = FALSE
    ) +
  geom_smooth(
    data = mu, aes(y = Estimate, ymin = Q2.5, ymax = Q97.5),
    stat = "identity", color = "black", alpha = 0.8, size = 1
    )

Modèle de régression, taille d’effet

Plusieurs méthodes pour calculer les tailles d’effet dans les modèles bayésiens. Gelman & Pardoe (2006) proposent une méthode pour calculer un \(R^{2}\) basé sur l’échantillon.

Modèle de régression, taille d’effet

\[ \begin{aligned} \rho^{2} &= \dfrac{\sum_{i = 1}^{n} \pi_{i}(\beta_{i} - \beta)^{2}}{\sigma^{2} + \sum_{i=1}^{n} \pi_{i}(\beta_{i} - \beta)^{2}} \\ \rho^{2} &= \dfrac{ \frac{1}{n} \sum_{i=1}^{n} \beta_{i}^{2}}{\sigma^{2} + \frac{1}{n} \sum_{i = 1}^{n} \beta_{i}^{2}} \\ \rho^{2} &= \dfrac{\beta^{2} \tau^{2}}{\sigma^{2} + \beta^{2} \tau^{2}}\\ \end{aligned} \]

post <- as_draws_df(x = mod4)
beta <- post$b_weight
sigma <- post$sigma

f1 <- beta^2 * var(d2$weight)
rho <- f1 / (f1 + sigma^2)

Modèle de régression, taille d’effet

posterior_plot(samples = rho, usemode = TRUE) + labs(x = expression(rho) )

summary(lm(height ~ weight, data = d2) )$r.squared

[1] 0.5696444

Modèle de régression, taille d’effet

bayes_R2(mod4)

   Estimate  Est.Error      Q2.5     Q97.5
R2 0.568369 0.02321423 0.5177365 0.6096139

bayes_R2(mod4, summary = FALSE)[, 1] %>%
    posterior_plot(usemode = TRUE) +
    labs(x = expression(rho) )

Résumé du cours

On a présenté un nouveau modèle à deux puis trois paramètres : le modèle gaussien, puis la régression linéaire gaussienne, permettant de mettre en relation deux variables continues.

Travaux pratiques - 1/2

Sélectionner toutes les lignes du jeu de données howell correspondant à des individus mineurs (age < 18). Cela devrait résulter en une dataframe de 192 lignes.

Travaux pratiques - 2/2

Imaginons que vous ayez consulté une collègue experte en allométrie (i.e., les phénomènes de croissance différentielle d’organes) et que cette dernière vous explique que ça ne fait aucun sens de modéliser la relation entre le poids et la taille… alors qu’on sait que c’est le logarithme du poids qui est relié à la taille !

Proposition de solution

# on garde seulement les individus ayant moins de 18 ans
d <- get(data(howell) ) %>% filter(age < 18)

priors <- c(
  prior(normal(150, 100), class = Intercept),
  prior(normal(0, 10), class = b),
  prior(exponential(0.01), class = sigma)
  )

mod7 <- brm(
  height ~ 1 + weight,
  prior = priors,
  family = gaussian(),
  data = d
  )

Proposition de solution

summary(mod7, prob = 0.89)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ 1 + weight 
   Data: d (Number of observations: 192) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-89% CI u-89% CI Rhat Bulk_ESS Tail_ESS
Intercept    58.27      1.42    56.02    60.53 1.00     3462     2903
weight        2.72      0.07     2.60     2.83 1.00     3359     2689

Family Specific Parameters: 
      Estimate Est.Error l-89% CI u-89% CI Rhat Bulk_ESS Tail_ESS
sigma     8.54      0.45     7.88     9.28 1.00     3724     3138

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Représenter les prédictions du modèle

# on crée un vecteur de valeurs possibles pour "weight"
weight.seq <- data.frame(weight = seq(from = 5, to = 45, length.out = 1e2) )

# on récupère les prédictions du modèle pour ces valeurs de poids
mu <- data.frame(
  fitted(mod7, newdata = weight.seq, probs = c(0.055, 0.945) )
  ) %>%
  bind_cols(weight.seq)

pred_height <- data.frame(
  predict(mod7, newdata = weight.seq, probs = c(0.055, 0.945) )
  ) %>%
  bind_cols(weight.seq)

# on affiche les 6 premières lignes de pred_height
head(pred_height)

  Estimate Est.Error     Q5.5    Q94.5   weight
1 71.80340  8.755395 58.00913 85.93774 5.000000
2 73.15025  8.619529 59.44546 86.84472 5.404040
3 74.00033  8.705761 60.03307 87.93907 5.808081
4 75.02742  8.526849 61.22442 88.70289 6.212121
5 76.36137  8.545242 62.87965 89.98825 6.616162
6 77.24049  8.606262 63.61413 91.13649 7.020202

Représenter les prédictions du modèle

d %>%
  ggplot(aes(x = weight, y = height) ) +
  geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
  geom_ribbon(
    data = pred_height, aes(x = weight, ymin = Q5.5, ymax = Q94.5),
    alpha = 0.2, inherit.aes = FALSE
    ) +
  geom_smooth(
    data = mu, aes(y = Estimate, ymin = Q5.5, ymax = Q94.5),
    stat = "identity", color = "black", alpha = 0.8, size = 1
    )

Proposition de solution

# on considère maintenant tous les individus
d <- howell

mod8 <- brm(
  # on prédit la taille par le logarithme du poids
  height ~ 1 + log(weight),
  prior = priors,
  family = gaussian(),
  data = d
  )

Proposition de solution

summary(mod8, prob = 0.89)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: height ~ 1 + log(weight) 
   Data: d (Number of observations: 544) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Population-Level Effects: 
          Estimate Est.Error l-89% CI u-89% CI Rhat Bulk_ESS Tail_ESS
Intercept   -23.59      1.35   -25.79   -21.43 1.00     3540     2776
logweight    47.02      0.38    46.41    47.64 1.00     3501     2956

Family Specific Parameters: 
      Estimate Est.Error l-89% CI u-89% CI Rhat Bulk_ESS Tail_ESS
sigma     5.16      0.16     4.91     5.42 1.00     3478     2686

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Représenter les prédictions du modèle

# on crée un vecteur de valeurs possibles pour "weight"
weight.seq <- data.frame(weight = seq(from = 5, to = 65, length.out = 1e2) )

# on récupère les prédictions du modèle pour ces valeurs de poids
mu <- data.frame(
  fitted(mod8, newdata = weight.seq, probs = c(0.055, 0.945) )
  ) %>%
  bind_cols(weight.seq)

pred_height <- data.frame(
  predict(mod8, newdata = weight.seq, probs = c(0.055, 0.945) )
  ) %>%
  bind_cols(weight.seq)

# on affiche les 6 premières lignes de pred_height
head(pred_height)

  Estimate Est.Error     Q5.5    Q94.5   weight
1 51.97455  5.233549 43.62893 60.35633 5.000000
2 57.36630  5.230244 48.90691 65.57163 5.606061
3 62.34145  5.178360 54.22963 70.82694 6.212121
4 66.66506  5.133290 58.60935 74.91199 6.818182
5 70.84443  5.184584 62.42015 78.99551 7.424242
6 74.29080  5.147859 66.11623 82.55444 8.030303

Représenter les prédictions du modèle

d %>%
  ggplot(aes(x = weight, y = height) ) +
  geom_point(colour = "white", fill = "black", pch = 21, size = 3, alpha = 0.8) +
  geom_ribbon(
    data = pred_height, aes(x = weight, ymin = Q5.5, ymax = Q94.5),
    alpha = 0.2, inherit.aes = FALSE
    ) +
  geom_smooth(
    data = mu, aes(y = Estimate, ymin = Q5.5, ymax = Q94.5),
    stat = "identity", color = "black", alpha = 0.8, size = 1
    )